Spatio-temporal models
Many real-world data sets vary in both space and time.
We often want to model how similarity (correlation) changes across both dimensions.
The good news 👍
The bad news 👎
fully spatio-temporal models are complex
model fitting can take a long time
different types of spatio-temporal data structures lead to different types of complexities
additional dependencies: in space and time
spatial and temporal behaviour can be independent on each other, or dependent: i.e. properties of the spatial structure vary through time (or vice versa)
Several options are available in inlabru
Areal Model
Geostatistical and Point Process models
We define a stochastic process:
\[ Z(s,t), \quad s \in \mathcal{S} \subset \mathbb{R}^d, \quad t \in \mathcal{T} \]
The covariance function:
\[ C((s_1,t_1),(s_2,t_2)) = \text{Cov}[Z(s_1,t_1), Z(s_2,t_2)] \]
Our goal: specify \(C(\cdot)\) to capture spatial, temporal, and joint dependencies.
Separable models \[ C((s_1,t_1),(s_2,t_2)) = C_S(s_1,s_2) \times C_T(t_1,t_2) \]
Interpretation:
Space and time act independently.
Spatial correlation not dependent on time lag.
Advantages: 😀
Simpler estimation.
Fewer parameters \(\rightarrow\) Easier computation.
Disadvantages: 😔
Non-Separable models
Separable models \[ C((s_1,t_1),(s_2,t_2)) = C_S(s_1,s_2) \times C_T(t_1,t_2) \]
Interpretation:
Space and time act independently.
Spatial correlation not dependent on time lag.
Advantages: 😀
Simpler estimation.
Fewer parameters \(\rightarrow\) Easier computation.
Disadvantages: 😔
Non-Separable models
\[ C((s_1,t_1),(s_2,t_2)) = f(|s_2-s_1|,|t_2-t_1|) \] Advantages: 😀
More realistic for dynamic physical processes.
Can model propagation or decay that varies over time.
Disadvantages: 😔
More parameters → complex estimation.
Need more data for estimation
Interpretation may be less straightforward.
We here focus on separable models
\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988 (simulated data!)
Space
\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988
Time
\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988
Time
The observation model
\[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]
We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)
Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)
Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.
The observation model \[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]
We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)
Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)
Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.
Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)
Here the spatial fields are different for each year, they are independent on each other.
The observation model
\[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]
We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)
Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)
Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.
Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)
Here the spatial fields are different for each year, they are independent on each other.
Model 3: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \[ u_{st} = \phi u_{st-1} + w_{st}, \text{ with } w_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T \]
Here we are adding more temporal dependence in the spatial term, using a AR1 type model Here the spatial fields are different for each year, they are independent on each other.
The observation model
\[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]
Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)
Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)
Model 3: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \[ u_{st} = \phi u_{st-1} + w_{st}, \text{ with } w_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T \]
Compare
Differences between Model 2 and 3 are not evident when estimating the past, but they would give different predictions
Model 3 in this case seem to better fit the data (more about these measures later)
model DIC WAIC
1 Model 1 14013.57 17381.01
2 Model 2 12346.99 12240.29
3 Model 3 11943.31 11923.44
mean sd 0.025quant 0.975quant
Precision for space 13.05 1.28 10.68 15.71
GroupRho for space 0.89 0.01 0.86 0.91
Precision for time 64.79 27.36 26.11 131.75
inlabruModels can be stuck together using either
replicate - the different slices are independent but share the hyperparametergroup - different dependence structures are implementedTypes of group model
[1] "exchangeable" "exchangeablepos" "ar1" "ar"
[5] "rw1" "rw2" "besag" "iid"
Note The group and replicate features can be used for more than space-time modeling!
(Simulated data)
The observation model
\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]
We are going to see 3 different models for the linear predictor \(\eta_{st}\)
Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)
Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.
Define mesh and spde model
border_simplified = st_simplify(border, dTolerance = 25)
mesh = fm_mesh_2d(boundary = border_simplified,
max.edge = c(30, 90), cutoff = 15,
offset = c(-0.05, -0.25),
crs = st_crs(border))
spde = inla.spde2.pcmatern(mesh = mesh,
prior.range = c(200, 0.5), # P(range < 100) = 0.5
prior.sigma = c(1, 0.5)) # P(sigma > 1) = 0.5Note
Coordinate system is in Km (this is recommended always!)
The size of the domain is ca \(630\times488 \text{ Km}^2\) we use this as a guideline to define the prior for the range.
Fit the model
cmp = ~ Intercept(1) +
time(time, model = "rw2") +
space(geometry, model = spde)
lik = bru_obs(formula = Y~ .,
data = dd)
fit_1 = bru(cmp, lik)
# plot time effect
p = fit_1$summary.random$time %>%
ggplot() +
geom_ribbon(aes(ID,
ymin = `0.025quant`,
ymax = `0.975quant`),
alpha = 0.5) +
geom_line(aes(ID, mean))Space-time Predictions
The observation model
\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]
We are going to see 3 different models for the linear predictor \(\eta_{st}\)
Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)
Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.
Model 2: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \(u_t(s)\sim\text{GF}(\rho,\sigma_u),\ t = 1,\dots,T\)
Here spatial fields are different for each year, are independent and share range and sd.
Fit the model
cmp = ~ Intercept(1) +
time(time, model = "rw2") +
space(geometry, model = spde,
replicate = time)
lik = bru_obs(formula = Y~ .,
data = dd)
fit_2 = bru(cmp, lik)
# plot time effect
p = fit_2$summary.random$time %>%
ggplot() +
geom_ribbon(aes(ID,
ymin = `0.025quant`,
ymax = `0.975quant`),
alpha = 0.5) +
geom_line(aes(ID, mean))Space-time Predictions
The observation model
\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]
We are going to see 3 different models for the linear predictor \(\eta_{st}\)
Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)
Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.
Model 2: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \(u_t(s)\sim\text{GF}(\rho,\sigma_u),\ t = 1,\dots,T\)
Here spatial fields are different for each year, are independent and share range and sd.
Model 3: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \[ u_t(s) = \phi\ u_{t-1}(s) + \omega(s),\text{ with } \omega(s)\sim\text{GF}(\rho,\sigma_u) \]
Here we are adding more temporal dependence in the spatial term, using a AR1 type model Here the spatial fields are different for each year, they are independent on each other.
Fit the model
h.spec <- list(rho = list(prior = 'pc.cor1',
param = c(0, 0.9)))
cmp = ~ Intercept(1) +
time(time, model = "rw2") +
space(geometry, model = spde,
group = time,
control.group = list(model = 'ar1',
hyper = h.spec))
lik = bru_obs(formula = Y~ .,
data = dd)
fit_3 = bru(cmp, lik)
# plot time effect
p = fit_3$summary.random$time %>%
ggplot() +
geom_ribbon(aes(ID,
ymin = `0.025quant`,
ymax = `0.975quant`),
alpha = 0.5) +
geom_line(aes(ID, mean))Space-time Predictions
Covariates that only vary in space
Covariates that vary in space and time
terra, SpatRaster object.class : SpatRaster
size : 147, 147, 1 (nrow, ncol, nlyr)
resolution : 4.465128, 3.11244 (x, y)
extent : -465.4573, 190.9165, 7031.885, 7489.413 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=23 +south +datum=WGS84 +units=km +no_defs
source(s) : memory
name : last
min value : -3.1660518
max value : 0.8164945
Setting up and running the model
Look at the results
Some covariates change in both space and time.
There are several ways of doing this….many quite confusing 🤪
Here is one recipe !
spatraster objectclass : SpatRaster
size : 147, 147, 12 (nrow, ncol, nlyr)
resolution : 4.465128, 3.11244 (x, y)
extent : -465.4573, 190.9165, 7031.885, 7489.413 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=23 +south +datum=WGS84 +units=km +no_defs
source(s) : memory
names : 1, 2, 3, 4, 5, 6, ...
min values : -3.1660518, -1.275877, -1.572957, -1.991793, -2.704403, -3.526074, ...
max values : 0.8164945, 2.442684, 2.221178, 2.395833, 4.012937, 1.116201, ...
Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -103.4191 ymin: 7460.922 xmax: -103.4191 ymax: 7460.922
Projected CRS: +proj=utm +zone=23 +south +datum=WGS84 +units=km +no_defs +type=crs
# A tibble: 3 × 3
Y time geometry
<dbl> <dbl> <POINT [km]>
1 1.19 1 (-103.4191 7460.922)
2 2.72 5 (-103.4191 7460.922)
3 2.79 6 (-103.4191 7460.922)
Note The layers names in the raster are the same as the elements of the column time in the data.
Implementation
mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept 0 0 0 0 1 0 0
cov -1 0 -1 -1 -1 -1 0
The .data. indicates the data that are input in the bru_obs() function
For space-time point processes one can also use the group and replicate features
The one difference is that one has to define the integration scheme in both space and time!
group feature)group feature)group feature)group feature)Where is our process defined (space and time)
Domain of integration
group feature)Results
inlabruinlabru. SeeLindgrenn et al., A diffusion-based spatio-temporal extension of Gaussian Matern > fields (2024) SORT 48